Thematische Karten erstellen mit R

Jan-Philipp Kolb

23 November 2017

Motivation - Deutschlands größte Klimasünder

Gliederung

Quellen für Polygonzüge

Hello World

library(maps)
map()

Das Paket maps - etwas detailierter

Grenzen sind recht grob:

map("world", "Germany")

Das Paket maps - Mehr Information

data(world.cities)
map("france")
map.cities(world.cities,col="blue")

Das Paket maptools

library(maptools)
data(wrld_simpl)
plot(wrld_simpl,col="royalblue")

Was sind shapefiles (.shp)?

head(wrld_simpl@data)
length(wrld_simpl)
## [1] 246
nrow(wrld_simpl@data)
## [1] 246

Einzelne Elemente des Datensatzes plotten

ind <- which(wrld_simpl$ISO3=="DEU")
plot(wrld_simpl[ind,])

wrld_simpl@data[ind,]
##     FIPS ISO2 ISO3  UN    NAME  AREA  POP2005 REGION SUBREGION   LON   LAT
## DEU   GM   DE  DEU 276 Germany 34895 82652369    150       155 9.851 51.11

Das R-Paket choroplethrMaps

library(ggplot2)
library(choroplethrMaps)
data(country.map)
ggplot(country.map, aes(long, lat, group=group)) + geom_polygon()

Eine Karte für die USA

data(state.map)
ggplot(state.map,aes(long,lat,group=group))+geom_polygon()

Andere Quellen für Shapefiles - Das Paket raster

library(raster)
LUX1 <- getData('GADM', country='LUX', level=1)
plot(LUX1)

Daten für das Luxemburg Beispiel

head(LUX1@data)
OBJECTID ID_0 ISO NAME_0 ID_1 NAME_1 HASC_1 CCN_1 CCA_1 TYPE_1 ENGTYPE_1 NL_NAME_1 VARNAME_1
1 131 LUX Luxembourg 1 Diekirch LU.DI NA District District Dikrech|Dikkrich
2 131 LUX Luxembourg 2 Grevenmacher LU.GR NA District District Gréivemaacher
3 131 LUX Luxembourg 3 Luxembourg LU.LU NA District District Lëtzebuerg|Luxemburg

Shapefiles bei Eurostat

BKG - Quelle für Kreise in Deutschland

library(maptools)
krs <- readShapePoly("vg250_ebenen/vg250_krs.shp")
plot(krs)
head(krs@data$RS)
## [1] 03401 03458 09473 05962 10046 05916
## 402 Levels: 01001 01002 01003 01004 01051 01053 01054 01055 01056 ... 16077

Die Kreise für Baden-Württemberg

BLA <- substr(krs@data$RS,1,2)
plot(krs[BLA=="08",])

Shapefiles für Wahlkreise

Ortsnetzbereiche

Quelle: Bundesnetzagentur

onb <- readShapePoly("onb_grenzen.shp")
head(onb@data)

Karte der Vorwahlbereiche

Einen größeren Vorwahlbereich ausschneiden

vwb <- as.character(onb@data$VORWAHL)
vwb1 <- substr(vwb, 1,2)
vwb7 <- onb[vwb1=="07",]
plot(vwb7)

vwb <- as.character(onb@data$ONB_NUMMER)
vwb1 <- substr(vwb, 1,1)
vwb7 <- onb[vwb1=="7",]
plot(vwb7)

Das Paket rgdal

library(rgdal)
## OGR data source with driver: ESRI Shapefile 
## Source: "post_pl.shp", layer: "post_pl"
## with 8270 features
## It has 3 fields
library(rgdal)
PLZ <- readOGR ("post_pl.shp","post_pl")

PLZ-Bereiche in Stuttgart

SG <- PLZ[PLZ@data$PLZORT99=="Stuttgart",]
plot(SG,col="chocolate1")

PLZ-Bereiche in Berlin

BE <- PLZ[PLZ@data$PLZORT99%in%c("Berlin-West","Berlin (östl. Stadtbezirke)"),]
plot(BE,col="chocolate2",border="lightgray")

Zwischenfazit - Quellen für Polygonzüge

Thematische Karten mit R erstellen

Thematische Karten erzeugen - das Paket sp

library(sp)
spplot(wrld_simpl,"POP2005")

Andere Einfärbungen wählen - Das Paket colorRamps

library(colorRamps)
spplot(wrld_simpl,"POP2005",col.regions=blue2red(100))

Es gibt auch noch mehr Farbverläufe im Paket colorRamps

spplot(wrld_simpl,"POP2005",col.regions=matlab.like(100))

Thematische Karten mit dem Paket choroplethr

library(choroplethr)
data(df_pop_state)
head(df_pop_state)

Eine thematische Karte mit choroplethr erstellen

state_choropleth(df_pop_state)

Nur drei Staaten darstellen

state_choropleth(df_pop_state,
                 title      = "2012 Population Estimates",
                 legend     = "Population",
                 num_colors = 1,
                 zoom       = c("california", "washington", "oregon"))

Eine Karte der US Counties

data(df_pop_county)
county_choropleth(df_pop_county)

Choroplethen Länder

data(df_pop_country)
country_choropleth(df_pop_country,
              title      = "2012 Population Estimates",
              legend     = "Population",
              num_colors = 1,
              zoom       = c("austria","germany",
                             "poland", "switzerland"))

Weltbank Daten

library(WDI) 
WDI_dat <- WDI(country="all", indicator=c("AG.AGR.TRAC.NO","TM.TAX.TCOM.BC.ZS"),
    start=1990, end=2000)
head(WDI_dat)

Weltkarte mit den Weltbank Daten

choroplethr_wdi(code="SP.DYN.LE00.IN", year=2012,
                title="2012 Life Expectancy")

Eurostat Daten

Sie können eine Statistik der Sparquote bei Eurostat downloaden.

http://ec.europa.eu/eurostat/web/euro-indicators/peeis

library(xlsx)
HHsr <- read.xlsx2("data/HHsavingRate.xls",1)

Zensus Ergebnisse

Zensus Atlas

https://ergebnisse.zensus2011.de/

Zensus Datenbank

Zensus Datenbank

Zensus Gemeindeergebnisse

zen <- read.csv2("data/Zensus_extract.csv")
# Personen mit eigener Migrationserfahrung
# mit beidseitigem Migrationshintergrund
zen2 <- data.frame(Personen_Mig=zen[,which(zen[9,]==128)],
                   Personen_Mig_bs=zen[,which(zen[9,]==133)])

datahub.io

Weltkulturerbestätten

url <- "https://raw.githubusercontent.com/Japhilko/
GeoData/master/2015/data/whcSites.csv"

whcSites <- read.csv(url) 
name_en date_inscribed longitude latitude area_hectares category states_name_fr
Cultural Landscape and Archaeological Remains of the Bamiyan Valley 2003 67.82525 34.84694 158.9265 Cultural Afghanistan
Minaret and Archaeological Remains of Jam 2002 64.51606 34.39656 70.0000 Cultural Afghanistan
Historic Centres of Berat and Gjirokastra 2005 20.13333 40.06944 58.9000 Cultural Albanie
Butrint 1992 20.02611 39.75111 NA Cultural Albanie
Al Qal’a of Beni Hammad 1980 4.78684 35.81844 150.0000 Cultural Algérie
M’Zab Valley 1982 3.68333 32.48333 665.0300 Cultural Algérie

Exkurs - OpenStreetMap Projekt

OpenStreetMap.org ist ein im Jahre 2004 gegründetes internationales Projekt mit dem Ziel, eine freie Weltkarte zu erschaffen. Dafür sammeln wir weltweit Daten über Straßen, Eisenbahnen, Flüsse, Wälder, Häuser und vieles mehr.

http://www.openstreetmap.de/

OSM Mapper

OSM Mapper

Export von OpenStreetMap Daten

<www.openstreetmap.org/export>

OpenStreetMap - Map Features

Overpass Turbo

https://overpass-turbo.eu/

https://overpass-turbo.eu/

Query Overpass

node
  [amenity=bar]
  ({{bbox}});
out;

Zwischenfazit - Quellen für inhaltliche Daten

Verknüpfung von Daten

Quelle: Geographic Information Systems and Remote Sensing

Daten verbinden - Beispiel Eurostat Daten

ind <- match(HHsr$geo,wrld_simpl@data$NAME)
ind <- ind[-which(is.na(ind))]
EUR <- wrld_simpl[ind,]
EUR@data$HHSR_2012Q3 <- as.numeric(as.character(HHsr[-(1:2),2]))
EUR@data$HHSR_2015Q2 <- as.numeric(as.character(HHsr[-(1:2),13]))

Karte mit Eurostat Indikator Household Saving Rate

spplot(EUR,c("HHSR_2012Q3","HHSR_2015Q2"))

Daten verbinden - Beispiel Bäckereien in Berlin

OSM als Datenquelle

(load("data/info_bar_Berlin.RData"))
## [1] "info"
addr.postcode addr.street name lat lon
79675952 13405 Scharnweberstraße Albert’s 52.56382 13.32885
86005430 NA NA Newton Bar 52.51293 13.39123
111644760 NA NA No Limit Shishabar 52.56556 13.32093
149607257 NA NA en passant 52.54420 13.41298
248651127 10115 Bergstraße Z-Bar 52.52953 13.39564
267780050 10405 Christburger Straße Immertreu 52.53637 13.42509

Verwendung des Pakets gosmd

devtools::install_github("Japhilko/gosmd")
library("gosmd")
pg_MA <- get_osm_nodes(object="leisure=playground","Mannheim")
pg_MA <- extract_osm_nodes(pg_MA,value='playground')

Matching

tab_plz <- table(info_be$addr.postcode)
ind <- match(BE@data$PLZ99_N,names(tab_plz))
ind
##   [1]  1  2  3  4  5  6  7  8 NA  9 NA NA NA NA NA 10 11 12 NA 13 14 15 16
##  [24] 17 18 19 20 21 22 23 24 25 NA 26 27 28 29 NA NA NA NA 30 NA 31 32 33
##  [47] 34 35 NA NA 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 NA 52 53
##  [70] NA 54 55 NA NA NA 56 57 58 59 60 NA NA NA NA NA 61 NA NA NA 62 NA NA
##  [93] NA NA NA NA NA NA NA 63 NA NA 64 NA 65 NA NA NA 66 NA NA NA NA 67 NA
## [116] NA 68 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [139] NA 69 70 NA 71 72 73 74 75 NA 76 NA NA NA NA NA NA NA NA NA NA NA NA
## [162] 77 NA 78 79 NA NA NA NA 80 NA NA NA NA 81 NA 82 83 84 NA NA NA NA NA
## [185] NA NA NA 85 NA NA

Daten anspielen

BE@data$num_plz <- tab_plz[ind]

Das Paket tmap

library(tmap)
BE@data$num_plz[is.na(BE@data$num_plz)] <- 0
qtm(BE,fill = "num_plz")

Mehr Informationen einbinden

load("data/osmsa_PLZ_14.RData")
PLZ99 PLZ99_N PLZORT99 nname EWZ_gem area_d EWZ_gemplz place_id osm_type osm_id lat lon display_name class type importance state city county plz2ort bakery bar biergarten butcher cafe chemist clothes college store food general cream kiosk mall pub restaurant supermarket population_density BLA gadmbla gadmkreis stop yes gadmgem gadmgemtyp gadmgem2 gadmgemtyp2 ort2plz ODdat zenEinw crossing bus_stop street_lamp traffic_signals land_cover.index land_cover.value land_cover.description elevation.value temp_Jan temp_Feb temp_Mar temp_Apr temp_May temp_Jun temp_Jul temp_Aug temp_Sep temp_Oct temp_Nov temp_Dez
0 01067 1067 Dresden Dresden, Stadt 512354 0.0008602 20494.16 144969068 relation 191645 51.0493286 13.7381437 Dresden, Sachsen, Deutschland place city 0.8162766 Sachsen Dresden Dresden 25 17 10 0 4 28 2 45 0 0 21 0 1 3 0 8 100 6 567 14 Sachsen Dresden 101 0 Dresden Einheitsgemeinde Dresden Stadt 1 0 4.0 121 48 162 87 22 Artificial surfaces and associated areas urban, water, vegetation, mountains, etc. 112 -0.7 0.4 3.9 8.4 13.3 16.9 18.5 18.0 14.3 9.8 4.4 1.0
1 01069 1069 Dresden Dresden, Stadt 512354 0.0006819 20494.16 144969068 relation 191645 51.0493286 13.7381437 Dresden, Sachsen, Deutschland place city 0.8162766 Sachsen Dresden Dresden 25 20 6 0 9 24 5 41 0 0 28 0 0 3 0 2 22 9 498 14 Sachsen Dresden 83 0 Dresden Einheitsgemeinde Dresden Stadt 1 0 5.0 113 40 105 96 22 Artificial surfaces and associated areas urban, water, vegetation, mountains, etc. 115 -0.8 0.3 3.8 8.4 13.3 16.7 18.4 17.9 14.4 9.9 4.4 1.0
2 01097 1097 Dresden Dresden, Stadt 512354 0.0004382 20494.16 144969068 relation 191645 51.0493286 13.7381437 Dresden, Sachsen, Deutschland place city 0.8162766 Sachsen Dresden Dresden 25 22 9 0 4 22 3 28 0 0 23 0 0 3 0 15 49 14 567 14 Sachsen Dresden 40 0 Dresden Einheitsgemeinde Dresden Stadt 1 0 5.5 98 20 33 55 22 Artificial surfaces and associated areas urban, water, vegetation, mountains, etc. 115 -0.7 0.3 3.8 8.4 13.3 16.7 18.4 18.0 14.4 9.9 4.5 1.0
3 01099 1099 Dresden Dresden, Stadt 512354 0.0067740 20494.16 144969068 relation 191645 51.0493286 13.7381437 Dresden, Sachsen, Deutschland place city 0.8162766 Sachsen Dresden Dresden 25 18 35 5 2 35 1 33 0 1 30 0 2 0 0 25 59 6 567 14 Sachsen Dresden 88 0 Dresden Einheitsgemeinde Dresden Stadt 1 0 0.0 38 41 24 37 4 Tree Cover, needle-leaved, evergreen urban, water, vegetation, mountains, etc. 250 -1.2 -0.3 3.1 7.6 12.5 16.0 17.6 17.4 13.7 9.3 3.8 0.4
4 01109 1109 Dresden Dresden, Stadt 512354 0.0034973 20494.16 144969068 relation 191645 51.0493286 13.7381437 Dresden, Sachsen, Deutschland place city 0.8162766 Sachsen Dresden Dresden 25 14 0 0 3 4 1 5 0 0 7 0 0 0 0 0 17 4 567 14 Sachsen Dresden 242 0 Dresden Einheitsgemeinde Dresden Stadt 1 0 1.5 47 119 230 58 22 Artificial surfaces and associated areas urban, water, vegetation, mountains, etc. 216 -1.0 -0.1 3.2 7.8 12.7 16.1 17.7 17.6 13.9 9.3 3.9 0.6
5 01127 1127 Dresden Dresden, Stadt 512354 0.0003626 20494.16 144969068 relation 191645 51.0493286 13.7381437 Dresden, Sachsen, Deutschland place city 0.8162766 Sachsen Dresden Dresden 25 6 1 0 3 4 0 6 0 0 6 0 0 0 0 5 13 3 567 14 Sachsen Dresden 44 0 Dresden Einheitsgemeinde Dresden Stadt 1 0 4.5 204 22 36 12 22 Artificial surfaces and associated areas urban, water, vegetation, mountains, etc. 112 -0.7 0.4 3.9 8.4 13.4 16.9 18.4 18.0 14.4 9.8 4.4 1.1

OSM-Daten - Bäckereien in Stuttgart

qtm(PLZ_SG,fill="bakery")

In welchem PLZ Bereich sind die meisten Bäckereien

kable(PLZ_SG@data[which.max(PLZ_SG$bakery),c("PLZ99","lat","lon","bakery")])
PLZ99 lat lon bakery
4964 70173 48.7784485 9.1800132 30

Das R-Paket RDSTK

library("RDSTK")

Die Daten für Stuttgart

PLZ_SG <- PLZ[PLZ@data$PLZORT99=="Stuttgart",]
Type_landcover Freq
Artificial surfaces and associated areas 26
Cultivated and managed areas 8
Tree Cover, needle-leaved, evergreen 1

Eine Karte der Flächenbedeckung erstellen

qtm(PLZ_SG,fill="land_cover.value")

Die Höhe in Stuttgart

qtm(PLZ_SG,fill="elevation.value")

Graphiken Stadtleben Stuttgart - das Paket ggmap

devtools::install_github("dkahle/ggmap")
install.packages("ggmap")

Eine erste Karte mit ggmap erzeugen

library(ggmap)
qmap("Stuttgart")

Karte für einen ganzen Staat

qmap("Germany")

Ein anderes zoom level

qmap("Germany", zoom = 6)

Karte für eine Sehenswürdigkeit

WIL <- qmap("Wilhelma",zoom=20, maptype="satellite")
WIL

ggmap - maptype satellite zoom 20

qmap('Stuttgart Hauptbahnhof', zoom = 15, maptype="hybrid")

Eine terrain Karte

qmap('Stuttgart Fernsehturm', zoom = 14,
 maptype="terrain")

ggmap - maptype watercolor

qmap('Stuttgart', zoom = 14,
 maptype="watercolor",source="stamen")

ggmap - source stamen

qmap('Stuttgart', zoom = 14,
 maptype="toner",source="stamen")

ggmap - maptype toner-lite

qmap('Stuttgart', zoom = 14,
 maptype="toner-lite",source="stamen")

# ggmap - maptype toner-hybrid

qmap('Stuttgart', zoom = 14,
 maptype="toner-hybrid",source="stamen")

ggmap - maptype terrain-lines

qmap('Stuttgart', zoom = 14,
 maptype="terrain-lines",source="stamen")

Geokodierung

Geocoding (…) uses a description of a location, most typically a postal address or place name, to find geographic coordinates from spatial reference data …

Wikipedia - Geocoding

library(ggmap)
geocode("Stuttgart")
##        lon      lat
## 1 9.182932 48.77585
lon lat
9.12036 48.81257

Reverse Geokodierung

Reverse geocoding is the process of back (reverse) coding of a point location (latitude, longitude) to a readable address or place name. This permits the identification of nearby street addresses, places, and/or areal subdivisions such as neighbourhoods, county, state, or country.

Quelle: Wikipedia

revgeocode(c(48,8))
## [1] "Unnamed Road, Somalia"

Die Distanz zwischen zwei Punkten

mapdist("Marienplatz Stuttgart","Hauptbahnhof Stuttgart")
##                    from                     to    m    km   miles seconds
## 1 Marienplatz Stuttgart Hauptbahnhof Stuttgart 3136 3.136 1.94871     488
##    minutes     hours
## 1 8.133333 0.1355556
mapdist("Marienplatz Stuttgart","Hauptbahnhof Stuttgart",mode="walking")
##                    from                     to    m    km    miles seconds
## 1 Marienplatz Stuttgart Hauptbahnhof Stuttgart 2505 2.505 1.556607    1874
##    minutes     hours
## 1 31.23333 0.5205556

Eine andere Distanz bekommen

mapdist("Marienplatz Stuttgart","Hauptbahnhof Stuttgart",mode="bicycling")
##                    from                     to    m    km    miles seconds
## 1 Marienplatz Stuttgart Hauptbahnhof Stuttgart 2722 2.722 1.691451     486
##   minutes hours
## 1     8.1 0.135

Geokodierung - verschiedene Punkte von Interesse

POI1 <- geocode("B2, 1 Mannheim",source="google")
POI2 <- geocode("Hbf Mannheim",source="google")
POI3 <- geocode("Mannheim, Friedrichsplatz",source="google")
ListPOI <-rbind(POI1,POI2,POI3)
POI1;POI2;POI3
##        lon      lat
## 1 8.462844 49.48569
##        lon      lat
## 1 8.469879 49.47972
##        lon      lat
## 1 8.475754 49.48304

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),
data = ListPOI)

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),col="red",
data = ListPOI)

ggmap - verschiedene Farben

ListPOI$color <- c("A","B","C")
MA_map +
geom_point(aes(x = lon, y = lat,col=color),
data = ListPOI)

ggmap - größere Punkte

ListPOI$size <- c(10,20,30)
MA_map +
geom_point(aes(x = lon, y = lat,col=color,size=size),
data = ListPOI)

Eine Route von Google maps bekommen

from <- "Mannheim Hbf"
to <- "Mannheim B2 , 1"
route_df <- route(from, to, structure = "route")

Mehr Information

http://rpackages.ianhowson.com/cran/ggmap/man/route.html

Eine Karte mit dieser Information zeichnen

qmap("Mannheim Hbf", zoom = 14) +
  geom_path(
    aes(x = lon, y = lat),  colour = "red", size = 1.5,
    data = route_df, lineend = "round"
  )

Wie fügt man Punkte hinzu

http://i.stack.imgur.com

Cheatsheet

https://www.rstudio.com/

Das Paket ggmap

library(ggmap)
lon_plz <- PLZ_SG@data[which.max(PLZ_SG$bakery),"lon"]
lat_plz <- PLZ_SG@data[which.max(PLZ_SG$bakery),"lat"]
mp_plz <- as.numeric(c(lon_plz,lat_plz))
qmap(location = mp_plz,zoom=15)

Das Paket osmar

library(osmar) 
src <- osmsource_api()
gc <- geocode("Stuttgart-Degerloch")
bb <- center_bbox(gc$lon, gc$lat, 800, 800)
ua <- get_osm(bb, source = src)
plot(ua)

Gebäude in diesem Ausschnitt plotten

bg_ids <- find(ua, way(tags(k=="building")))
bg_ids <- find_down(ua, way(bg_ids))
bg <- subset(ua, ids = bg_ids)
bg_poly <- as_sp(bg, "polygons")  
plot(bg_poly)

Gebäude und Straßen im Ausschnitt

Schlussfolie

Vielen Dank für die Aufmerksamkeit

Erweiterungsmöglichkeiten

http://rpubs.com/Japhilko82/Rleaflet

Karte auf Carto zur Energieerzeugung